import kotlinx.cinterop.*
import platform.posix.*
import kotlin.system.getTimeNanos
import kotlin.math.round

fun main(args: Array<String>) {
    // 运行24次基准测试
    for (i in 1..24) {
        val startTime = getTimeNanos()
        main_(args)
        val endTime = getTimeNanos()

        // 计算并打印耗时
        val duration = (endTime - startTime).toDouble() / 1_000_000_000.0
        println("duration: ${formatDouble(duration, 3)}")
    }
}

fun main_(args: Array<String>) {
    var n = 100
    if (args.isNotEmpty()) {
        n = args[0].toIntOrNull() ?: 100
    }

    val result = SampleSpectralNorm().approximate(n)
    println(formatDouble(result, 9))
}

// 修复浮点数格式化函数，使用循环计算幂次
fun formatDouble(value: Double, decimalPlaces: Int): String {
    // 用循环计算10的幂次，替代pow函数
    var multiplier = 1.0
    repeat(decimalPlaces) {
        multiplier *= 10.0
    }

    // 明确类型转换，解决重载歧义
    val rounded = round(value * multiplier) / multiplier

    return if (decimalPlaces == 0) {
        rounded.toInt().toString()
    } else {
        val parts = rounded.toString().split('.')
        val integerPart = parts[0]
        val decimalPart = if (parts.size > 1) parts[1] else "0"
        if (decimalPart.length >= decimalPlaces) {
            "$integerPart.${decimalPart.substring(0, decimalPlaces)}"
        } else {
            "$integerPart.${decimalPart.padEnd(decimalPlaces, '0')}"
        }
    }
}

class SampleSpectralNorm {
    fun approximate(n: Int): Double {
        // 创建单位向量
        val u = DoubleArray(n) { 1.0 }
        val v = DoubleArray(n) { 0.0 }

        // 幂法迭代20步
        for (i in 0 until 10) {
            multiplyAtAv(n, u, v)
            multiplyAtAv(n, v, u)
        }

        // 计算 v.Bv/(v.v) 的平方根，其中 B = A^T A
        var vBv = 0.0
        var vv = 0.0
        for (i in 0 until n) {
            vBv += u[i] * v[i]
            vv += v[i] * v[i]
        }

        return kotlin.math.sqrt(vBv / vv)
    }

    /** 计算无限矩阵A的元素(i,j) */
    private fun a(i: Int, j: Int): Double {
        val idx = (i + j) * (i + j + 1) / 2 + i + 1
        return 1.0 / idx
    }

    /** 向量v乘以矩阵A */
    private fun multiplyAv(n: Int, v: DoubleArray, av: DoubleArray) {
        for (i in 0 until n) {
            av[i] = 0.0
            for (j in 0 until n) {
                av[i] += a(i, j) * v[j]
            }
        }
    }

    /** 向量v乘以矩阵A的转置 */
    private fun multiplyAtv(n: Int, v: DoubleArray, atv: DoubleArray) {
        for (i in 0 until n) {
            atv[i] = 0.0
            for (j in 0 until n) {
                atv[i] += a(j, i) * v[j]
            }
        }
    }

    /** 向量v乘以矩阵A，再乘以矩阵A的转置 */
    private fun multiplyAtAv(n: Int, v: DoubleArray, atAv: DoubleArray) {
        val u = DoubleArray(n)
        multiplyAv(n, v, u)
        multiplyAtv(n, u, atAv)
    }
}
